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1  Introduction 

The  study  of  detonation  waves  dates  back  to  the  late  19th  century,  where 
Chapman  [1]  and  Jouguet  [2]  modeled  detonations  as  a  shock  wave  supported 
by  the  heat  release  of  the  combustible  material  in  an  infinitely  thin  zone, 
where  all  chemistry  and  diffusive  transport  takes  place.  Later  Zel’dovich  [3], 
von  Neumann  [4],  and  Boring  [5]  independently  represented  the  detonation 
as  a  confluence  of  a  shock  wave  moving  at  a  detonation  velocity,  followed 
by  a  chemical  reaction  zone  of  finite  length;  this  came  to  be  known  as  the 
ZND  model  for  a  detonation  wave. 

While  the  true  structure  of  detonation  waves  inevitably  calls  for  multi¬ 
dimensional  effects,  the  simple  ID  structure  still  provides  a  rich  spectrum 
of  dynamical  features  which  are  worthy  of  detailed  exploration.  This  is 
especially  important  for  the  study  of  deflagration  to  detonation  transitions 
[6]  and  sustained  oscillating  or  galloping  detonations  [7,  8] .  For  a  spark- 
induced  detonation,  as  the  detonation  decays  towards  the  self-sustaining 
Chapman-Jouguet  mode  from  an  over-driven  mode,  one  obtains  a  sequence 
of  physical  oscillations  between  the  flame  and  shock  front.  The  numerical 
analysis  of  this  effect  has  been  explored  by  Cambier  [9]  using  highly-resolved 
numerical  simulations,  albeit  with  only  a  2"'Lorder  shock  capturing  scheme. 
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Similar  calculations  of  the  dynamics  of  detonation  oscillations,  albeit  with 
a  reduced  chemistry  model,  have  been  performed  elsewhere  [10-12]. 

In  the  present  study,  we  utilize  high-order  spatially  accurate  numerical 
methods  in  order  to  achieve  grid  convergence  and  to  reduce  or  eliminate 
numerical  diffusion  effects  while  providing  a  detailed  analysis  of  the  non¬ 
linear  dynamics  involved  in  resolving  detonations  with  complex  reaction 
kinetics.  The  dynamical  characteristics  and  coupling  of  large  and  small 
length  scale  physics  associated  with  detonation  instabilities  is  explored  in 
detail. 


2  Numerical  Methodology 

2.1  Governing  Equation 

Since  we  are  interested  in  the  dynamics  associated  with  the  inviscid  problem, 
we  solve  the  Euler  equations  for  a  multi-species,  real  gas  with  a  chemical 
source  term,  as  shown  below: 

kQ  +  |.F(Q)  =  s(Q)  (1) 

where  the  vectors  represented  by  Q,  F,  and  S  are,  respectively: 


Q=  ,F  = 

pu 

\e 


pu^  -|-  P 
{E  +  P)u 


The  total  mixture  density  \s  p  =  ps-,  with  ps  the  density  of  species  s,  and 
P  and  u  represent  pressure  and  velocity,  respectively.  The  total  energy  E 
may  be  written  as 


E  =  y  pCv{T)dT  +  (3) 

where  c„(r)  =  Yhs  PsCvs{T)/p  is  the  mass-averaged  specific  heat  at  constant 
volume.  Note  that  the  total  energy,  Eq.  (3),  does  not  include  the  potential 
chemical  energy  of  the  mixture,  given  by  Y2s  Ps^Os,  where  cqs  is  the  formation 
energy  of  the  species.  Therefore,  the  formation  energy  appears  in  Eq. 
(2)  as  a  source  term,  to  account  for  the  exo-  or  endo-thermicity  of  the 
chemical  reactions  and  guarantee  energy  conservation.  This  approach  has 
been  found  preferable  to  including  the  formation  energy  in  the  definition  of 
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the  mixture’s  total  energy  -  in  which  case  no  heat  of  formation  term  is  added 
to  the  total  energy,  due  to  chemical  reactions  is  necessary  -  as  it  can  yield 
better  accuracy  at  composition  discontinuities.  In  the  governing  equation, 

Eq.  (2),  ijjs  represents  the  net  production  rate  of  the  species: 

(hs  =  ^  Vrskfr  p'j"'  -  ^  Vrshr  p'i' 
r  j  r  j 

Urs  —  Vj-g  —  Vyg 

where  and  are  the  coefficients  of  species  in  the  forward  and 
backward  reactions,  respectively,  and  kfr  and  k^r  and  the  forward  and  back¬ 
ward  chemical  rates  of  the  reaction. 

It  is  noted  that  some  recent  studies  (eg.  Romick  et  al  [13])  suggest 
that  the  inclusion  of  physical  diffusion  in  a  detonation  simulation  is  needed 
to  reduce  or  eliminate  the  effects  of  numerical  diffusion  on  resolution  re¬ 
quirements  and  to  generate  more  physically  accurate  instabilities.  Without 
these  physical  dissipation  mechanisms  of  diffusion  or  viscosity,  the  chemistry 
provides  the  only  time  and  length  scales,  and  the  spectrum  of  dynamical 
features  of  the  system  can  potentially  include  very  high-frequency  phenom¬ 
ena,  especially  if  the  chemical  kinetics  are  very  stiff.  The  extent  to  which 
these  phenomena  can  be  reproduced  is  at  question  here,  since  an  attempt  at 
reproducing  the  physical  characteristics  of  the  problem  would  require  multi¬ 
dimensional  calculations.  Therefore,  it  is  essential  that  we  eliminate  as  much 
as  possible  the  artificial  scales  of  numerical  dissipation  in  order  to  capture 
the  chemically-induced  dynamics.  For  this  purpose,  in  addition  to  the  use 
of  high-order,  low-dissipation  numerical  methods,  we  have  placed  a  special 
emphasis  on  grid  resolution  until  grid-independent  results  are  obtained. 

2.2  Numerical  Scheme  and  Chemical  Reaction  Mechanism 

In  the  present  study,  a  5*^-order  accurate  Monotonicity  Preserving  Scheme(MP5) 
as  formulated  by  Suresh  and  Huynh  [14]  is  principally  used  for  high  order 
interpolation  of  a  system  of  governing  equations.  The  MP  schemes  form 
a  class  of  high-order  accurate,  shock-capturing  numerical  algorithms  espe¬ 
cially  equipped  for  solving  non-linear  hyperbolic  systems  of  conservation  law 
equations.  The  MP5  scheme  has  a  5*^  order  accuracy  in  smooth  regions  of 
flow- field,  while  reducing  to  1^*  order  at  flow  discontinuities,  a  necessary  con¬ 
dition  to  avoid  spurious  oscillations.  This  makes  MP  schemes  well  suited 
for  resolving  flow- fields  where  shocks  and  flame  fronts  are  present.  The 
MP  schemes  interpolate  the  value  of  the  conserved  variables  (listed  in  the 
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definiton  of  the  Q  vector)  at  the  cell  interface.  The  interface  flux  is  then 
solved  utilizing  the  Roe  Flux  (RF)  [15]  building  block  with  the  Harten,  Lax, 
van  Leer,  and  Einfeldt  [16]  entropy  fix  applied  to  the  eigenvalues.  The  RF 
scheme  is  slightly  more  computationally  expensive  than  one  of  its  popular 
variants.  Local  Lax  Friedrich  (LLF),  but  is  less  dissipative,  making  it  opti¬ 
mal  for  studying  stability  and  resolving  high-frequency  wave  structures.  A 
S^Aorder  Total  Variation  Diminishing  (TVD)  Runge-Kutta  time  integration 
scheme  is  used  in  conjunction  with  MP5. 

A  fifth  order  spatially  accurate  Weighted  Essentially  Non-Oscillatory 
(WENO),  scheme  as  formulated  by  Shu  &  Osher  [17]  was  used  in  congunc- 
tion  with  the  Advection-Diffusion-Reaction  (ADER)  scheme  of  Titerev  & 
Toro  [18]  to  form  the  spatially  5*^  order  accurate  ADER-WEN05  or  sim¬ 
ply  AW5  scheme.  ADER  schemes  utilize  the  high  order  spatial  derivatives 
calculated  by  the  underlying  scheme,  e.g.  WENO,  to  generate  the  tempo¬ 
ral  derivatives  using  Rosonov’s  procedure  [19],  for  example,  dtu  =  \dxU^ 
dttu  =  XdxxU,  etc.  With  the  high  order  temporal  derivatives,  a  simple  Tay¬ 
lor  series  expansion  is  performed  to  acquire  a  higher  order  temporally  and 
spatially  accurate  scheme.  Since  the  expensive  Runge-Kutta  time  integra¬ 
tion  steps  are  no  longer  required,  ADER- WENO  is  extremely  efficient  and 
well  suited  for  parallel  computation. 

Since  complex  chemical  reaction  mechanisms  are  simulated,  an  operator 
splitting  [20]  approach  is  implemented  to  facilitate  the  efficient  computa¬ 
tion  of  the  coupling  with  the  chemical  kinetics.  The  latter  generally  form 
a  stiff  system  and  for  stability  reasons,  are  computed  with  a  point-implicit 
backwards  Euler  method.  As  in  a  previous  study  [9],  the  chemical  kinetics 
of  a  diluted  hydrogen-oxygen-nitrogen  mixture  were  solved.  The  chemistry 
includes  eight  reacting  species,  H2,  O2,  H,  O,  OH,  HO2,  H2O2,  H2O,  and 
the  non-reacting  diluent  V2.  Thirty  eight  elementary  reactions  are  used  in 
this  mechanism  and  the  backward  rates  are  computed  from  equilibrium  con¬ 
stants.  Eor  this  chemical  system  of  moderate  size,  the  Gaussian  elimination 
scheme  is  the  optimal  choice  for  solving  the  chemical  kinetics  implicitly. 

It  should  be  emphasized  that  the  operator-splitting  approach  used  here  is 
not  the  only  choice  available  to  us.  In  fact,  a  more  accurate  coupling  would 
be  obtained  by  solving  a  fully-coupled  system  (convection  and  chemistry 
together)  with  an  implicit  Runge-Kutta  method.  However,  the  computa¬ 
tional  cost  would  increase  at  least  3-fold  on  an  already  challenging  problem, 
even  in  ID.  Eurthermore,  we  are  mostly  concerned  with  the  reproduction  of 
wave  phenomena  within  a  small  distance,  i.e.  the  induction  length;  our  grid 
spacing  is  sufficiently  small  that  the  time  step  becomes  dominated  by  the 
Courant-Eriedrichs-Lewy  (CEL)  limitation  rather  than  the  significant  evo- 
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lution  of  the  chemical  kinetics.  Thus,  we  leave  the  comparison  of  coupling 
and  time-stepping  methods  as  a  topic  for  future  work. 

2.3  Code  Validation  and  Testing 

In  order  to  ensure  that  small  scale  features  are  properly  captured  and  the 
leading  shock  of  the  detonation  is  properly  resolved  by  the  numerical  scheme, 
a  series  of  validation  test  were  performed.  The  first  case,  as  prescribed  by 
Shu  &  Osher  in  [21]  demonstrates  how  well  a  numerical  scheme  resolves  a 
shock  propagating  through  a  entropy  disturbance.  Figure  1  shows  an  ex¬ 
ample  of  the  spatial  variation  in  density  at  a  given  instant  of  time  for  this 
problem.  It  is  clear  that  the  MP5  scheme  is  able  to  capture  a  range  of 
different  complex  wave  structures  in  the  flow  quite  accurately.  The  second 
case,  referred  to  as  the  blast- wave  problem  in  [22],  tests  the  ability  of  a 
numerical  scheme  to  capture  the  interaction  of  strong  shocks  with  shocks 
and  expansion  waves.  Figure  2  similarly  demonstrates  that  complex  waves 
interactions  are  accurately  captured  with  modest  grid  resolution.  By  com¬ 
parison,  the  AW5  solution  shown  in  3  is  slightly  more  diffusive.  This  is 
despite  the  fact  that  no  Runge-Kutta  time  stepping  is  used  in  that  case. 
Since  both  methods  are  formally  5*^-order  accurate  in  smooth  regions,  the 
difference  in  quality  of  the  solutions  is  due  to  the  limiter.  In  the  following, 
most  of  the  results  regarding  the  detonation  dynamics  will  be  obtained  with 
the  MP5  scheme,  while  the  AW5  scheme  is  used  to  evaluate  the  influence 
of  numerical  accuracy.  We  should  also  emphasize  that  the  results  shown 
in  Figures  1-3  do  not  use  any  artificial  compression  at  the  contact  discon¬ 
tinuities  (CD).  Thus,  the  profile  near  x  ~  0.600  in  the  blast-wave  problem 
could  be  made  much  sharper;  however,  the  non-linear  procedure  necessary 
for  sharpening  CDs  may  also  artificially  sharpen  the  gradients  in  acoustic 
and  entropy  waves  and  affect  the  dynamics  within  the  induction  region  of  a 
detonation;  to  guarantee  the  absence  of  numerical  artifacts  in  the  solution, 
no  artificial  compression  is  used  in  the  computations. 

In  the  present  study,  we  have  modeled  the  detailed  kinetics  of  the  sim¬ 
plest  combustion  system.  Another  approach,  commonly  chosen  in  funda¬ 
mental  studies  of  detonation  dynamics,  is  a  1-step  model,  in  which  the  entire 
chemistry  is  described  by  the  evolution  of  a  single  progress  variable  that  fol¬ 
lows  an  exponential  relaxation  with  a  characteristic  time-scale  given  by  the 
induction  delay.  This  progress  variable  is  also  associated  with  the  fractional 
amount  of  heat  released  into  the  flow.  In  that  model,  the  delay  follows  a 
simple  exponential  fit  r*  ~  The  delay  being  essentially  caused  by 

the  need  for  a  sufficient  amount  of  radicals  from  chain-branching  reactions. 
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and  the  production  of  those  being  an  endo-thermic  process,  the  parameter 
Ea  in  this  formulation  is  an  averaged  activation  energy  of  the  key  radical- 
producing  reactions.  This  is  a  reasonable  approximation  to  the  chemistry 
in  that  region,  albeit  within  limits.  We  have  used  the  detailed  chemistry 
to  compute  and  parametrize  the  induction  delay  as  a  function  of  initial 
temperature,  and  pressure  (the  mixture  being  held  fixed  to  stoichiometric 
hydrogen-air).  As  shown  in  Figure  4,  the  delay  does  follow  an  exponential 
form  ti  oc  as  expected.  The  parameter  (3  in  this  formulation  is 

an  averaged  activation  energy  of  the  key  radical-producing  reactions.  How¬ 
ever,  this  approach  yields  unrealistic  profiles  of  the  post-shock  region,  since 
the  heat  release  is  gradual.  A  better  description  is  obtained  with  the  2-step 
model  [12],  where  the  heat  release  is  associated  with  a  second  progress  vari¬ 
able,  whose  evolution  can  start  only  at  the  end  of  the  induction  delay,  which 
now  follows  a  linear  time  variation.  While  this  allows  a  separation  between 
the  induction  and  heat  release  zones,  this  model  fails  in  several  ways.  First, 
the  heat  release  is  assumed  independent  of  temperature,  which  is  unrealistic, 
as  the  flow  heating  accelerates  the  combustion;  generally  speaking,  a  stiff 
differential  equation  for  the  progress  variable  can  be  used  to  reproduce  this 
non-linear  effect,  but  the  dynamics  can  be  different  from  the  real  conditions. 
Second,  when  the  flame  is  accelerated  towards  the  shock,  the  two  reaction 
zones  (induction  and  flame)  start  to  merge,  even  if  species  diffusion  is  ne¬ 
glected;  the  enforcing  of  two  separate  zones  with  a  2-step  model  modifies  the 
dynamics  of  the  strongly  coupled  shock- flame  system.  Thus,  we  have  used  a 
realistic  chemical  system  in  order  to  gain  a  better  understanding  of  the  true 
dynamics;  the  comparison  with  the  results  obtained  with  a  simplified  2-step 
model  would  then  allow  us  to  correlate  the  results  with  specific  features  of 
the  chemical  kinetics. 

3  Detonation  Dynamical  Phenomena 

3.1  Ignition  and  Instabilities 

Direct  initiation  of  the  ID  detonation  in  the  study  is  obtained  in  a  chamber 
filled  with  a  stoichiometric  mixture  of  H2  and  air  (temperature  300  K  and 
pressure  1  atm)  by  setting  a  region  adjacent  to  an  end-wall  of  the  simu¬ 
lated  shock  tube  at  high  pressure  (40  atm)  and  temperature  (1500  K),  as 
a  simulated  spark.  The  direction  initiation  is  preferable  to  a  deflagration- 
to  detonation  transition  (DDT),  since  the  latter  is  much  more  sensitive  to 
initial  conditions  and  grid  resolution,  requires  species  diffusion,  and  a  very 
long  computational  domain  to  reproduce  both  the  DDT  and  the  subsequent 
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evolution  of  the  detonation.  Nevertheless,  even  direct  initiation  is  sensitive 
to  initial  conditions  and  resolution.  The  requirements  to  achieve  detonation 
ignition  with  the  MP5  scheme  include  a  grid  cell  size  of  less  than  Ax  =  50 
//m  and  a  distributed  simulated  spark  region  (of  a  length  0.25  to  0.5  cm) 
with  sufficiently  high  presure.  Figure  5(a)  illustrates  the  pressure  contours 
of  a  spark  ignited  mixture  which  does  not  achieve  detonation,  while  Fig¬ 
ure  5(b)  illustrates  contours  of  the  same  mixture  and  grid  resolution  with 
a  higher  spark  pressure  that  achieves  detonation.  The  initial  spark  con¬ 
ditions  can  be  related  to  the  minimum  energy  and  kernel  size  for  direct 
initiation  [23].  Similar  studies  by  He  &  Karagozian  [24]  indicate  there  are 
both  pressure  and  temperature  requirements  for  the  spark.  Once  a  satisfac¬ 
tory  kernel/spark  size  and  pressure  were  found,  these  initial  conditions  were 
kept  for  the  remainder  of  the  studies. 

The  succesful  detonation  initiation  event  proceeds  in  two  phases;  first, 
the  gas  in  the  spark  region  rapidly  burns  and  increases  the  pressure  to 
even  larger  values,  in  a  nearly  constant-volume  combustion  process.  This 
high  pressure  generates  a  strong  shock  which  propagates  into  the  unburnt 
mixture,  which  itself  is  ignited  after  a  time  delay  and  rapidly  burns,  starting 
from  the  region  closest  to  the  spark.  In  a  scenario  described  as  the  SWACER 
mechanism  [28],  the  combustion  wave  overtakes  the  leading  shock  and  the 
coalescence  of  the  two  fronts  leads  to  extremely  high  peak  pressures  for  a 
very  short  time.  This  event  is  easily  identified  in  the  trace  of  the  peak 
pressure  versus  time,  as  shown  in  Figure  6,  and  is  referred  to  hereafter  as 
the  “re-explosion”  event  (the  first  explosion  taking  place  within  the  initial 
spark  region).  Two  different  grid  sizes  are  used  in  Figures  6(a)  &  6(b).  A 
more  detailed  examination  for  the  dynamics  for  these  two  different  grid  sizes 
are  shown  in  Figures  7  &  8,  as  will  be  discussed  below. 

The  high  pressure  of  this  second  explosion  event  initiates  another  strong 
shock,  followed  after  an  an  induction  length  {£,  in  the  reference  frame  of  the 
shock)  by  the  combustion  zone.  The  flame  is  initially  strongly  coupled  to 
the  shock  {£  0)  and  the  wave  is  over-driven,  i.e.,  its  speed  exceeds  that 

of  the  Chapman-Jouget  (CJ)  detonation.  As  the  degree  of  overdrive  decays 
and  the  detonation  approaches  the  CJ  limit,  instabilities  begin  to  appear, 
as  shown  in  6(a)  after  35  ^s  and  6(b)  after  25  fis. 

Because  the  re-explosion  is  clearly  identified,  we  can  use  this  feature  to 
conduct  a  preliminary  study  of  the  effect  of  grid  resolution.  Thus,  we  exam¬ 
ined  the  variation  of  the  measured  time  delay  to  the  second  explosion  texp, 
for  given  initial  spark  conditions.  In  this  study,  the  uniform  grid  spacing 
Ax  was  varied  from  l^m  to  20/um.  Figure  9  illustrates  how  the  time  to 
re-explosion  for  both  the  MP5  and  AW5  schemes  varies  with  the  grid  res- 
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olution.  The  peak  in  texp  for  the  MP5  scheme  curve  in  Figure  9  occurs  at 
approximately  7.5^m,  whereas  it  is  slightly  lower,  at  approximately  5iJ,m, 
for  the  AW5  scheme.  Since  the  latter  is  more  diffusive  than  MP5,  this  is 
not  surprising  that  the  curve  exhibits  a  smoother  profile.  The  most  striking 
feature  here  is  the  non-monotonic  behavior,  i.e.  the  presence  of  a  maxi¬ 
mum  time  to  explosition,  which  delineates  two  regimes.  For  grid  resolutions 
Ax  below  this  critical  value,  the  simulation  is  in  a  “convectively”  dominant 
regime  where  the  combination  of  the  numerical  scheme  and  grid  resolution 
is  sufficient  to  mitigate  the  effects  of  numerical  diffusion  to  the  detonation 
formation.  Both  schemes  produce  similar  values  of  texp  for  Ax  <  2^m, 
but  for  this  complex  kinetics  scheme,  time  convergence  for  texp  may  not  be 
reached  except  for  Ax  <  l^um,  consistent  with  Powers  &  Paolucci’s  find¬ 
ings  [25].  For  grid  resolutions  greater  than  the  critical  value  in  Figure  9, 
we  enter  the  numerically  dissipative  regime  where  coupling  of  the  fluid  and 
kinetics  is  enhanced  due  to  numerical  diffusion  of  temperature  and  chemical 
concentrations. 

Given  this  result,  we  conclude  that  the  spacing  should  be  less  than  5 
or  7.5  fim,  depending  on  the  numerical  scheme.  To  obtain  results  truly 
independent  of  the  numerical  effects,  we  should  have  Ax  ^  0,  since  the  two 
methods  converge  -  presumably,  since  we  are  of  course  unable  to  run  such  a 
calculation  ~  towards  a  single  value  for  texp  in  that  limit.  If  we  allow  ourselves 
a  given  error  level  (say  5%),  we  can  use  a  spacing  of  approximately  5  ;um. 
In  these  studies  of  the  detonation  propagation  and  shock-flame  coupling 
dynamics,  most  of  the  results  were  obtained  for  a  spacing  of  Ax  =  2.5/um, 
thus  suggesting  good  accuracy.  Nevertheless,  the  time  to  re-explosion  is  but 
one  parameter  that  is  affected  by  grid  resolution,  and  calculations  of  the 
shock-flame  instabilities  were  repeated  for  several  grid  spacings  to  confirm 
the  accuracy  of  the  results. 

The  present  spark-ignited  detonation  simulations,  utilizing  the  two  high- 
order  schemes  mentioned  above,  demonstrate  the  appearance  of  different 
instability  modes.  Figure  6  illustrates  the  peak  pressure  of  in  the  entire 
computational  domain,  for  two  different  grid  sizes.  After  the  re-explosion, 
the  detonation  is  strongly  overdriven,  as  mentioned  earlier;  this  is  charac¬ 
terized  as  a  rather  smooth  region  in  the  peak  pressure  trace,  between  10 
and  25'*“^s,  in  the  high-resolution  case  of  Figure  6(b).  Instabilities  appear 
when  the  detonation  becomes  close  to  the  CJ  condition,  starting  with  a 
small- amplitude,  but  high-frequency  mode  -  hereafter  referred  to  as  “HF” 
For  the  low-resolution  case,  the  transition  to  the  high-frequency  (HF)  mode 
occurs  at  a  time  of  the  order  of  30 //s,, shown  in  detail  in  Figure  7(a).  For 
the  high-resolution  case,  the  transition  occurs  earlier,  at  25  fis.  Around  45 


/is,  ther  is  a  gradual  shift  towards  a  lower  frequency  but  high- amplitude 
mode  -  referred  to  here  as  “HA”.  This  transition  is  gradual  in  the  high- 
resolution  case,  with  both  modes  coexisting  during  a  period  of  time  (betwen 
approximately  44  and  48  ^s),  as  shown  in  more  detail  in  Figure  8(b).  In  the 
low-resolution  case  of  8(a),  however,  the  behavior  is  less  gradual  and  more 
chaotic.  Clearly,  there  are  significant  effects  of  the  grid  resolution  on  the 
dynamics  of  the  instabilities;  furthermore,  the  appearance  of  sharp  features 
in  the  traces  also  suggests  that  special  care  must  be  exercised  in  avoiding 
numerical  procedures  which  can  arbitrarily  sharpen  gradients,  as  mentioned 
above.  Furthermore,  there  are  clear  differences  in  the  specific  dynamical 
features  of  these  instabilities,  e.g.,  in  the  time  of  initiation  of  high  frequency 
behavior  and  in  the  temporal  waveforms  of  instabilities. 

The  nature  of  these  features  bears  more  quantitative  exploration.  A  sim¬ 
ple  fast  Fourier  transform  (FFT)  can  be  used  to  find  the  spectral  content  of 
these  two  modes  and  verify  their  grid  independence,  as  shown  in  Figure  10. 
The  spectral  content  of  the  simulations  is  relatively  insensitive  to  the  grid 
spacing,  as  long  as  we  are  below  the  critical  value  for  the  start  of  numerical 
diffusion.  For  both  HA  instabilities  near  a  frequency  of  approximately  0.35 
MHz  and  HF  instabilities  near  2.3  MHz,  there  is  relatively  little  difference 
in  results  for  Ax  <  2.5/xm.  At  very  high  frequencies,  above  4  MHz,  there 
is  some  grid  dependency.  Frequencies  above  4  MHz  are  not  seen  at  coarse 
resolutions,  and  the  spectral  content  around  4.5  MHz  is  likely  to  be  a  har¬ 
monic  of  the  strong  2.3  MHz  signal.  Note  that  since  the  instabilities  develop 
for  a  finite  time  only,  the  sampling  statistics  of  the  FFT  are  limited.  Using 
wavelet  decomposition  did  not  provide  improvements  in  the  signal-to-noise 
ratio. 

These  results  confirm  earlier  findings  [9]  as  well  as  more  recent  studies 
[10,11],  which  have  shown  that  a  detonation  at  the  CJ  limit  has  two  distinct 
instability  modes.  The  high  frequency  mode  always  appears  first  and  marks 
the  transition  from  a  ‘stable’  CJ  detonation.  These  fluctuations  in  key 
properties  (e.g.  species  concentration,  temperature,  and  pressure)  of  the 
fluid  within  the  induction  zone  are  described  by  Oran  and  Boris  [26]  as  ‘hot 
spots’.  In  the  present  study,  we  have  found  that  these  ‘hot  spots’  contribute 
to  an  initial  stage  of  the  flame  dynamics.  In  this  regime,  the  induction  length 
is  very  small  {£  «  icj):  and  acoustic  waves  generated  by  the  perturbed 
chemistry  are  rapidly  transmitted  to  the  shock.  Because  there  is  a  very 
limited  amount  of  fluid  that  can  participate  in  the  fluctuation  of  the  heat 
release,  only  low-amplitude  perturbations  of  the  CJ  peak  pressure  appear. 
As  these  acoustic  waves  reach  the  leading  shock  and  strengthen  it,  their 
frequency  can  be  measured  as  that  of  the  fluctuations  of  the  peak  pressure. 
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Eventually  the  average  induction  length  continues  to  increase  and  a  sec¬ 
ond  mode  can  be  seen,  which  directly  couples  the  flame  speed  with  the 
shock,  resulting  in  fluctuations  with  lower  frequency  but  much  higher  am¬ 
plitude,  as  will  be  discussed  below.  While  similar  studies  by  Ng  et  al  [27], 
observed  stability  limits  of  detonations  using  complex  kinetics,  were  were 
able  to  further  resolve  the  detonation  with  higher  schemes  and  more  points. 


3.2  Simplified  Model 

A  model  for  the  induction  zone  can  be  constructed,  as  composed  of  a  leading 
shock,  a  heated,  post-shock  medium(fluid),  and  a  flame  front,  all  of  which 
are  illustrated  in  Figure  11.  A  single  period  of  the  detonation  oscillation 
can  be  described,  in  the  rest-frame  of  the  shock,  as  follows.  Fluctuations 
at  the  flame  front  create  an  acoustic  (pressure)  disturbance,  which  travels 
at  the  acoustic  wave  speed,  Xac,  through  the  induction  zone  until  it  reaches 
the  leading  shock.  Upon  contact,  the  pressure  fluctuation  carried  by  the 
acoustic  wave  will  accelerate  the  shock  and  therefore  alter  the  post-shock 
conditions,  thus  creating  an  entropy  disturbance  (temperature  fluctuation). 
This  entropy  disturbance  will  propagate  back  into  the  induction  zone  at  the 
entropy  wave  speed,  Xen,  toward  the  flame  front.  Upon  contact  with  the 
flame,  the  entropy  wave  will  create  a  new  acoustic  disturbance  in  the  flame, 
and  the  cycle  will  repeat.  Figure  11  illustrates  this  phenomenon,  while  the 
equations  below  quantify  the  entropy  and  acoustic  wave  speeds. 


Xenix j 


U2{x,t) 


(5) 


Xac{x,  t) 


c{x,t)  -  U2{x,t) 


(6) 


where  u{x,  t)  is  the  fluid  velocity,  c{x,  t)  is  local  speed  of  sound,  D{t)  is  the 
detonation  velocity,  and  U2{x,t)  =  \u{x,t)  —  D{t)\  is  the  post-shock  fluid 
velocity  in  the  detonation  reference  frame. 

From  the  wave  speeds,  the  period  of  the  cycle,  r,  can  be  expressed  by 


r  = 


Xac{x,t) 


dx  + 


-^f 


Xenix,  t) 


dx 


(7) 


where  Xg  =  (t  —  to)  ■  D{t)  is  the  position  of  the  shock  and  xj  is  the  position 
of  the  flame. 

At  a  0*^-order  approximation,  the  fluid  properties  in  the  induction  re¬ 
gion,  Q2{x,t),  are  assumed  to  weakly  vary  with  time  for  a  given  half  cycle. 
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~  0  &  —  0.  From  this  approximation,  the  period  can  be 

determined 


Ca^b  “1“  ^a^b  ^6— >c  ^b—>c 

where  i  is  the  period  averaged  induction  length,  a  ^  b  is  the  fluid  state  at 
the  acoustic  wave  half-cycle,  b  ^  c  is  the  fluid  state  at  the  entropy  wave 
half  cycle,  and  Mq,  c^,  and  are  the  the  fluid  speed  in  the  detonation 
reference  frame,  speed  of  sound,  and  average  detonation  speed  for  half-cycle 
a,  respectively.  The  model  for  acoustic  and  entropy  half  cycles  are  illustrated 
in  Figure  12,  in  correspondence  to  observed  oscillations.  From  the  period  of 
the  combined  cycles,  the  frequency  in  oscillations  of  the  peak  pressure  trace 
is  /  =  T-^ 

3.3  Discussion  &;  Analysis 

From  the  simplified  model  expressed  in  Eq.  (8),  data  were  extracted  at  dif¬ 
ferent  peak  pressure  cycles  from  the  high  frequency  as  well  as  high  amplitude 
regime.  Figure  13  illustrates  the  evolution  of  the  induction  zone  tempera¬ 
ture  profile  in  the  detonation  reference  frame  for  a  given  period  of  the  high 
amplitude  mode.  Using  the  simulation  induction  zone  data  and  the  period 
from  Eq.  (8),  the  frequency  /  ~  310  khz  is  estimated,  which  is  in  excel¬ 
lent  agreement  with  that  obtained  from  the  spectral  analysis  (310  ±  40  fe/iz) 
shown  in  Eigure  10.  Performing  the  same  analysis  on  the  high  frequency 
mode,  illustrated  in  Eigure  14,  the  frequency  /  2.08  MHz  is  estimated, 

which  is  in  good  agreement  with  the  spectral  analysis  (2.29±0.4  Mhz)  which 
is  shown  in  Eigure  10. 

The  ‘hot  spot’,  which  appears  in  the  high  frequency  mode  temperature 
profiles,  Eigure  14,  is  of  particular  interest.  It  appears  to  inhibit  the  flames 
progress  toward  the  flame  front  by  pre-burning  the  fluid  in  the  induction 
zone.  This  pre-ignition  effect  does  not  allow  the  flame  to  accelerate,  and 
in  the  case  of  the  high  frequency  mode  the  detonation  is  still  strongly  over¬ 
driven.  The  combination  of  the  short  induction  length,  due  to  the  strong 
overdrive,  and  the  inhibited  flame  lead  to  the  high  frequency  and  relatively 
low  amplitude  oscillations  in  the  peak  pressure  trace. 

In  the  high  amplitude  regime,  Eigure  13,  the  flame  is  unihibited  by 
the  presence  of  the  ‘hot  spot’,  which  allows  it  to  propagate  through  the 
induction  zone  unmolested.  A  SWACER-like  [28]  mechanism  appears  to 
govern  the  high  amplitude  regime  where  in  the  flame  is  furthest  from  the 
leading  shock,  imax,  it  is  uninhibited  by  the  ‘hot  spot’  and  will  burn  the  fluid 
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in  the  induction  and  accelerate  toward  the  shock,  releasing  large  amounts 
energy. 


4  Conclusion  Sz.  Remarks 

The  present  studies  indicate  that  it  is  possible  to  utilize  a  spatially  high  or¬ 
der  accurate  scheme  with  appropriate  resolution  of  induction  (approximately 
100  points  per  induction  half  length)  to  explore  preferred  re-explosion  and 
instability  modes  of  a  ID  finite  rate  complex  kinetics  detonation.  Appropri¬ 
ate  resolution  of  the  unstable  detonation  with  complex  kinetics  requires  an 
examination  of  a  number  of  parameters  in  addition  to  peak  pressure,  includ¬ 
ing  dominant  frequencies,  time  to  re-explosion,  and  surely  other  phenomena 
as  well.  A  simple  model  for  the  transmission  of  acoustic  and  entropy  waves 
creating  the  different  instability  modes  (high  frequency,  low  amplitude  or 
HF  and  low  frequency,  high  amplitude  or  HA)  replicates  very  well  the  dom¬ 
inant  modes  and  appears  to  be  a  reasonable  description  for  the  mechanism 
by  which  the  instabilities  arise. 

There  are  a  number  of  open  questions  that  remain  in  the  exploration 
of  unstable  detonations  with  complex  reaction  kinetics.  A  comparison  of 
results  from  the  present  studies  with  those  using  a  simplified  (e.g.,  two- 
step)  model  would  help  to  ascertain  the  importance  of  the  chemistry  in 
such  overdriven  detonations.  Comparison  of  results  with  reduced  chemistry 
models  for  a  lower  overdrive  is  similarly  important.  The  question  of  the 
required  degree  of  grid  resolultion  is  also  relevant  to  these  issues.  These  and 
other  phenomena  will  be  explored  in  future  studies. 
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Figure  1:  Computed  (symbols)  and  exact  solutio  of  density  profile  in  Shu- 
Osher  test  problem  at  a  given  time.  The  MP5  scheme  is  used  here. 
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Figure  2:  Computed  (symbols)  and  exact  solution  of  density  profile  in  Blast- 
Wave  test  problem  at  a  given  time.  The  MP5  scheme  is  used  here. 
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Figure  3:  Computed  (symbols)  and  exact  solution  of  density  profile  in  Blast- 
Wave  test  problem  at  a  given  time.  The  AW5  scheme  is  used  here. 


Figure  4:  Induction  Delay  versus  initial  temperature  for  various  pressures. 

0{P) 

Curves  can  be  fitted  as  t  ~  a{P)e  t 
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(a)  Pspark  =  20atm  with  0.25cm  spark  length  where  detonation  is  not  achieved 
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(b)  Pspark  =  50atm  with  0.5cm  spark  length  where  detonation  is  achieved 

Figure  5:  Pressure  contour  of  a  spark  ignited  H2-A1T  mixture  with  Ax  = 
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(a)  Ax  =  12.5  fim 


(b)  Ax  =  2.5  pm 

Figure  6:  Peak  pressure  time  history  of  a  spark-ignited  ff2-air  mixture  sim¬ 
ulated  with  two  different  grid  cell  sizes  Ax.  The  time  to  re-explosion,  texp, 
high  amplitude  mode,  HA,  and  high  frequency  mode,  HF,  are  illustrated. 
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(a)  Ax  =  12.5  fim 


Figure  7:  High  frequency  portion  of  peak  pressure  time  history  as  in  Fig.  6, 
simulated  with  two  different  grid  cell  sizes  Ax. 
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Figure  8:  High  amplitude  portion  of  peak  pressure  time  history  as  in  Fig. 
6,  simulated  with  two  different  grid  cell  sizes  Ax. 
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Ax  [fim] 


Figure  9:  Time  to  re-explosion  versus  grid  resolution  for  both  MP5  (with 
3rd  Qj-fjgj-  Runge-Kutta)  and  AW5  schemes. 


Figure  10:  Fourier  transform  of  peak  pressure  traces  for  the  unstable  deto¬ 
nation  with  3  different  sizes  of  Ax. 
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Figure  11:  Model  of  shock,  flame,  induction  zone,  and  transmission  of  en¬ 
tropy  and  acoustic  waves  for  shock-flame  coupling. 


Figure  12:  Simplified  peak  pressure  cycle(left)  used  to  represent  the  numer¬ 
ically  observed  pulsations  in  peak  pressure(right). 
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Figure  13:  Temperature  distribution  in  shock  reference  frame  at  different 
times  within  the  acoustic  wave  (a)  and  entropy  wave  (b)  cycle.  Data  is 
extracted  from  high  amplitude  mode  results. 
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Figure  14:  Temperature  distribution  in  shock  reference  frame  at  different 
times  within  the  entropy  wave  (a)  and  acoustic  wave  (b)  cycle.  Data  is 
extracted  from  high  frequency  mode  results. 
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